A Graph-Theoretical Approach for 
the Analysis and Model Reduction of 
Complex-Balanced Chemical Reaction Networks 

Shodhan Rao* Arjan van der Schaft^ Bayu Jayawardhana"'" 

November 29, 2012 

Abstract 

In this paper we derive a compact mathematical formulation describing the dynamics of chemical 
reaction networks that are complex-balanced and are governed by mass action kinetics. The formulation 
is based on the graph of (substrate and product) complexes and the stoichiometric information of these 
complexes, and crucially uses a balanced weighted Laplacian matrix. It is shown that this formulation 
leads to elegant methods for characterizing the space of all equilibria for complex-balanced networks and 
for deriving stability properties of such networks. We propose a method for model reduction of complex- 
balanced networks, which is similar to the Kron reduction method for electrical networks and involves the 
computation of Schur complements of the balanced weighted Laplacian matrix. 

Keywords: Weighted Laplacian matrix, linkage classes, zero-deficiency networks, persistence conjecture, 
equilibria, Schur complement. 

1 Introduction 

The analysis and control of large-scale chemical reaction networks poses a main challenge to bio-engineering 
and systems biology. Chemical reaction networks involving several hundreds of chemical species and reactions 
are common in living cells [22]. The complexity of their dynamics is further increased by the fact that the 
chemical reaction rates are intrinsically nonlinear. In particular mass action kinetics, which is the most 
basic form of expressing the reaction rates, corresponds to differential equations which are polynomial in 
the concentrations of the chemical species. Chemical reaction network dynamics thus are a prime example 
of complex networked dynamical systems, and there is a clear need to develop a mathematical framework 
for handling their complexity. 
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One of the issues in formalizing complex chemical reaction network dynamics is the fact that their 
graph representation is not immediate; due to the fact that chemical reactions usually involve more than 
one substrate chemical species and more than one product chemical species. This problem is resolved by 
associating the complexes of the reactions, i.e. the left-hand (substrate) and right-hand (product) sides of 
each reaction, with the vertices of a graph, and the reactions with the edge3^} The resulting directed graph, 
called the graph of complexes, is characterized by its incidence matrix B. Furthermore the stoichiometric 
matrix S of the chemical reaction network, expressing the basic balance laws of the reactions, can be 
factorized as S = ZB, with the complex- stoichiometric matrix Z encoding the expressions of the complexes 
in the various chemical species. We note that an alternative graph formulation of chemical reaction networks 
is the species-reaction graph E] H] , which is a bipartite graph with one part of the vertices corresponding 
to the species and the remaining part to the reactions, and the edges expressing the involvement of the species 
in the reactions. 

Using the graph of complexes formalism, we have developed in [29] [23] a compact mathematical formula- 
tion for a class of mass-action kinetics chemical reaction networks, which is characterized by the assumption 
of the existence of an equilibrium for the reaction rates; a so-called thermodynamic equilibrium. This 
corresponds to the thermodynamically justified assumption of microscopic reversibility, with the resulting 
conditions on the parameters of the mass action kinetics usually referred to as the Wegscheider conditions. 
The resulting class of mass action kinetics reaction networks are called (detailed) -balanced. A main feature 
of the formulation of |29j is the fact that the dynamics of a detailed-balanced chemical reaction network 
is completely specified by a symmetric weighted Laplacian matrix, defined by the graph of complexes and 
the equilibrium constants, together with an energy function, which is subsequently used for the stability 
analysis of the network. In particular, the resulting dynamics is shown [24J to bear close similarity with 
consensus algorithms for symmetric multi-agent systems. (In fact, it is shown in |29j that the so-called 
complex- affinities asymptotically reach consensus.) Furthermore, as shown in |17j . the framework can be 
readily extended from mass action kinetics to (reversible) Michaelis-Menten reaction rates. 

On the other hand, the assumption of existence of a thermodynamical equilibrium requires reversibility 
of all the reactions of the network, while there are quite a few well-known irreversible chemical reaction 
network models, including the McKeithan network to be explained shortly afterwards. Motivated by such 
examples we will extend in this paper the results of |29] by considering the substantially larger class of 
complex-balanced reaction networks. A chemical reaction network is called complex-balanced if there exists a 
vector of species concentrations at which the combined rate of outgoing reactions from any complex is equal 
to the combined rate of incoming reactions to the complex, or in other words each of the complexes involved 
in the network is at equilibrium. The notion of complex-balanced networks was first introduced in |16| and 
studied in detail in |1 1|. IT5l [7J [25| [9] . These systems have also been called as toric dynamical systems in the 
literature (see [7]). 

An example of a complex-balanced network is the model of T-cell interactions due to |19j (see also |27j ) 

lr This approach is originating in the work of Horn and Jackson [16]; see also Othmer |21| and [3] for nice exposes and 
additional insights. 



depicted in Figure [T] This chemical reaction network model arises in immunology and was proposed by 
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Figure 1: McKeithan's network 

McKeithan in order to explain the selectivity of T-cell interactions. With reference to Figure [TJ T and 
M represent a T-cell receptor and a peptide-major histocompatibility complex (MHC) respectively and 
T + M is a complex for the network. For i = 1, . . . , N, Ci represent various intermediate complexes in 
the phosphorylation and other intermediate modifications of the T-cell receptor T; k p ^ represents the rate 
constant of the i th step of the phosphorylation and k-±^ is the dissociation rate of the i th complex. In the 
following, we denote by [A] the concentration of a species A participating in a chemical reaction network. 
The governing law of the reaction network is the law of mass action kinetics. This leads to the following 
set of differential equations describing the rate of change of concentrations of various species involved in the 
network: 
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Observe that if the left hand side of each of the above equations is set to zero, all the concentrations [Cj] for 
i = 0, . . . , N—l, can be parametrized in terms of [Cn] . Since all the rate constants and dissociation constants 
are positive, it is easy to see that there exists a set of positive concentrations {[T], [M], [Co], . . . , [Cat]} for 
which the right hand sides of the equations ([T]) vanish. This implies that McKeithan's network is complex- 
balanced. On the other hand, all reactions in this network are irreversible, and thus McKeithan's network 
is not detailed-balanced. 

The main aim of this paper is to show how the compact mathematical formulation of detailed-balanced 
chemical reaction networks derived in |29| can be extended to complex-balanced networks (such as McKei- 
than's network). Indeed, the crucial difference between detailed-balanced and complex-balanced networks 
will turn out to be that in the latter case the Laplacian matrix is not symmetric anymore, but still bal- 
anced (in the sense of the terminology used in graph theory and multi-agent dynamics). In particular, 
complex-balanced chemical networks will be shown to bear close resemblance with asymmetric consensus 
dynamics with balanced Laplacian matrix. Exploiting this formulation it will be shown how the dynamics 



of complex-balanced networks share important common characteristics with those of detailed-balanced net- 
works, including a similar characterization of the set of all equilibria and the same stability result stating 
that the system converges to an equilibrium point uniquely determined by the initial condition. For the 
particular case when the complex-stoichiometric matrix Z is injective as in the McKeithan's network case 



(see Remark 4.7), the same asymptotic stability results were obtained before in |27j . Furthermore, while 
for detailed-balanced networks it has been shown in (29] that all equilibria are in fact thermodynamic equi- 
libria in this paper the similar result will be proved that all equilibria of a complex-balanced network are 
complex-equilibria. Similar results have already been proved in [16]; however, the proofs presented in the 
current paper are much more concise and insightful as compared to those presented in [16] . 

Furthermore, based on our formulation of complex-balanced networks exhibiting a balanced weighted 
Laplacian matrix associated to the graph of complexes, we will propose a technique for model-reduction of 
complex-balanced networks. This technique is similar to the Kron reduction method for model reduction of 
resistive electrical networks described in [15] : see also [TU1 [28] . Our technique works by deleting complexes 
from the graph of complexes associated with the network. In other words, our reduced network has fewer 
complexes and usually fewer reactions as compared to the original network, and yet the behavior of a number 
of significant metabolites in the reduced network is approximately the same as in the original network. Thus 
our model reduction method is useful from a computational point of view, specially when we need to deal 
with models of large-scale chemical reaction networks. Mathematically our approach is based on the result 
that the Schur complement (with respect to the deleted complexes) of the balanced weighted Laplacian 
matrix of the full graph of complexes is again a balanced weighted Laplacian matrix corresponding to the 
reduced graph of complexes. 

The paper is organized as follows. In Section 2, we introduce tools from stoichiometry of reactions and 
graph theory that are required to derive our mathematical formulation. In Section 3, we explain mass action 
kinetics, define complex-equilibria and complex-balanced networks and then derive our formulation. In 
Section 4, we derive equilibrium and stability properties of complex-balanced networks using our formulation. 
In Section 5, we propose a model reduction method for complex-balanced networks, while Section 6 presents 
conclusions based on our results. 

Notation: The space of n dimensional real vectors is denoted by R n , and the space of m x n real matrices 
by R mxn . The space of n dimensional real vectors consisting of all strictly positive entries is denoted by W} 
and the space of n dimensional real vectors consisting of all nonnegative entries is denoted by W] . The rank 
of a real matrix A is denoted by rank A dim(V) denotes the dimension of a set V. Given a±, . . . , a n £ R, 
diag(ai, . . . ,a n ) denotes the diagonal matrix with diagonal entries oi, . . . ,a n ; this notation is extended to 
the block-diagonal case when a\,...,a n are real square matrices. Furthermore, ker A and span A denote the 
kernel and span respectively of a real matrix A. If U denotes a linear subspace of R m , then U 1 - denotes 
its orthogonal subspace (with respect to the standard Euclidian inner product). l m denotes a vector of 
dimension m with all entries equal to 1. The time-derivative j£(t) of a vector x depending on time t will be 
usually denoted by x. 



Define the mapping Ln : M™ — > M m , x h-> Ln(x), as the mapping whose i-th component is given as 
(Ln(x))j := ln(xj). Similarly, define the mapping Exp : R' m — > M™, x t— >■ Exp(x), as the mapping whose 
i-th component is given as (Exp(x)) i := exp(xj). Also, define for any vectors x,z € M. m the vector x ■ z € M m 
as the element-wise product (x • z\ := XiZi, i = 1,2, ... ,m, and the vector | 6 W 71 as the element-wise 
quotient (|)^ := jr, i = l,--- ,m. Note that with these notations Exp(x + z) = Exp(x) • Exp(z) and 
Ln(x • z) = Ln(x) + Ln(z), Ln (f ) = Ln(x) — Ln(z). 

2 Chemical reaction network structure 

In this section, we introduce the tools necessary in order to derive our mathematical formulation of the 
dynamics of complex-balanced networks. First we introduce the concept of stoichiometric matrix of a 
reaction network. We then define the concept of a graph of complexes, which was first introduced in the 
work of Horn k, Jackson and Feinberg ( [HI HSJ [T2J, [13] ) . 

2.1 Stoichiometry 

Consider a chemical reaction network involving m chemical species (metabolites), among which r chemical 
reactions take place. The basic structure underlying the dynamics of the vector x £ W± of concentrations 
Xi, i = 1, • • • , m, of the chemical species is given by the balance laws x = Sv, where S is an m x r matrix, 
called the stoichiometric matrix. The elements of the vector v S W are commonly called the (reaction) 
fluxes or rates. The stoichiometric matrix S, which consists of (positive and negative) integer elements, 
captures the basic conservation laws of the reactions. It already contains useful information about the 
network dynamics, independent of the precise form of the reaction rate v(x). Note that the reaction rate 
depends on the governing law prescribing the dynamics of the reaction network. 

We now show how to construct the stoichiometric matrix for a reaction network with the help of an 
example shown in Figure [2} 

X A 

X 1 + 2X 2 - X 3 2X 1 + X 2 
Figure 2: Example of a reaction network 

Note that the above reaction has 3 irreversible and 1 reversible reaction leading to in total 5 reactions 
among four species X\, X 2 , X3 and X4. Since the stoichiometric matrix maps the space of reactions to 
the space of species, it has dimension 4x5. The entry of S corresponding to the i th row and j th column 
is obtained by subtracting the number of moles of the i th species on the product side from that on the 



substrate side for the j reaction. Thus 
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for the reaction network depicted in Figure [2j 

In this paper, we focus only on closed chemical reaction networks meaning those without external fluxes. 
Therefore unless otherwise mentioned, our chemical reaction networks do not have any external fluxes. 

If there exists an m-dimensional row-vector k such that kS = 0, then the quantity kx is a conserved 
quantity or a conserved moiety for the dynamics x = Sv(x) for all possible reaction rates v = v(x). Indeed, 
kSv(x) 



dt KX 



0. Later on, in Remark 



3.3 



it will be shown that law of conservation of mass leads to a 
conserved moiety of a chemical reaction network. 

Note that for all possible fluxes the solutions of the differential equations x = Sv{x) starting from an 
initial state xq will always remain within the affine space 

S XQ := {x G | x- x G imS}. (2) 

S XQ has been referred to as the positive stoichiometric compatibility class (corresponding to xq) in p~3 | [25]. [Tj . 



2.2 Graph of complexes 

The structure of a chemical reaction network cannot be directly captured by an ordinary graph. Instead, 
we will follow an approach originating in the work of Horn and Jackson [16], introducing the space of 
complexes. The set of complexes of a chemical reaction network is simply defined as the union of all the 
different left- and righthand sides (substrates and products) of the reactions in the network. Thus, the 
complexes corresponding to the network ^ are X\ + 2X2, X3, 2X\ + X2 and X4. 

The expression of the complexes in terms of the chemical species is formalized by an m x c matrix Z, 
whose a-th column captures the expression of the a-th complex in the m chemical species. For example, for 
the network depicted in Figure [2j 

10 2 
2 10 

Z = 

10 
1 

We will call Z the complex- stoichiometric matrix of the network. Note that by definition all elements of the 
matrix Z are non- negative integers. 

Since the complexes are the left- and righthand sides of the reactions, they can be naturally associated 
with the vertices of a directed graph Q with edges corresponding to the reactions. Formally, the reaction 
a — > f3 between the a-th (reactant) and the /3-th (product) complex defines a directed edge with tail vertex 



being the a-th complex and head vertex being the /3-th complex. The resulting graph will be called the 
graph of complexes. 

Recall, see e.g. [3], that any graph is defined by its incidence matrix B. This is a c x r matrix, c being 
the number of vertices and r being the number of edges, with (a,j)-th element equal to —1 if vertex a is 
the tail vertex of edge j and 1 if vertex a is the head vertex of edge j, while otherwise. 

Obviously, there is a close relation between the matrix Z and the stoichiometric matrix S. In fact, it is 
easily checked that 

S = ZB, hence x = ZBv(x) (3) 
with B the incidence matrix of the graph of complexes. 



3 The dynamics of complex-balanced networks governed by mass action 
kinetics 

In this section, we first recall the dynamics of species concentrations of reactions governed by mass action 
kinetics. We then define complex-balanced networks and derive a compact mathematical formulation for 
their dynamics. 

3.1 The general form of mass action kinetics 

Recall that for a chemical reaction network, the relation between the reaction rates and species concentrations 
depends on the governing laws of the reactions involved in the network. In this section, we explain this 
relation for reaction networks governed by mass action kinetics. In other words, if v denotes the vector of 
reaction rates and x denotes the species concentration vector, we show how to construct v (x). The reaction 
rate of the j-th reaction of a mass action chemical reaction network, from a substrate complex Sj to a 
product complex Vj, is given as 

m z 

1=1 

where Zi p is the (i, p)-th element of the complex-stoichiometric matrix Z, and kj > is the rate constant 
of the j-th. reaction. Without loss of generality we will assume throughout that for every j, the constant kj 
is positive. 

This can be rewritten in the following way. Let Zs i denote the column of the complex-stoichiometric 
matrix Z corresponding to the substrate complex of the j-th reaction. Using the mapping Ln : W\ — > R c 
as defined at the end of the Introduction, the mass action reaction equation Q for the j-th reaction from 
substrate complex Sj to product complex Vj can be rewritten as 

Vj(x) = kj exp (ZjTn(x)). (5) 



Based on the formulation in ([5]), we can describe the complete reaction network dynamics as follows. Let 
the mass action rate for the complete set of reactions be given by the vector v(x) = v\(x) ■ ■ ■ v r (x) 
For every a, tt G {1, . . . , c}, define 

C* a := {j e{l,...,r}\ (a,n) = (S j ,V j )} 

and a ncr := X^'eC^ ^j- Thus if there is no reaction a — > ir, then a n(7 = 0. Define the weighted adjacency 
matrix A of the graph of complexes as the matrix with (tt, cr)-th element a na , where it, a £ {l,...,c}. 
Furthermore, define L := A — A, where A is the diagonal matrix whose (p,p)-th element is equal to the 
sum of the elements of the p-th column of A. Let B denote the incidence matrix of the graph of complexes 
associated with the network. By definition of L, we have \ T C L = 0. It can be verified that the vector 
Bv(x) for the mass action reaction rate vector v(x) is equal to — LExp (Z T Ln(x)^ , where the mapping 
Exp : M c — > Ml has been defined at the end of the Introduction. Hence the dynamics can be compactly 
written as 

x = -ZLExp (Z T Ln(x)) (6) 

A similar expression of the dynamics corresponding to mass action kinetics, in less explicit form, was already 
obtained in [27] , 

3.2 Complex-balanced networks 

We now define a class of reaction networks known as complex-balanced networks. This class was first defined 
in the work of Horn and Jackson (see p. 92 of [H]). We first define a complex-equilibrium of a reaction 
network. 

Definition 3.1. Consider a chemical reaction network with dynamics given by the equation Q. A vector 
of concentrations x* £ is called a complex-equilibrium if Bv(x*) = 0. Furthermore, a chemical reaction 
network is called complex-balanced if it admits a complex- equilibrium. 

It is easy to see that any complex-equilibrium is an equilibrium for the network, but the other way 
round need not be true (since Z need not be injective). We now explain the physical interpretation of a 
complex-equilibrium. Observe that 

Bv{x*) = (7) 

consists of c equations where c denotes the number of complexes. Among all the reactions that the i th 
complex Cj is involved in, let Oi denote the set of all the reactions for which Cj is the substrate complex 
and let 2j denote the set of all reactions for which C{ is the product complex. The i th of equations ^ can 
now be written as 



It follows that at a complex-equilibrium, the combined rate of outgoing reactions from any complex is equal 
to the combined rate of incoming reactions to the complex. In other words, at a complex-equilibrium, every 
complex involved in the network is at equilibrium. 

In |15j and [9], conditions have been derived for a chemical reaction network governed by mass action 
kinetics to be complex-balanced. 

Remark 3.2. A thermodynamically balanced or detailed-balanced chemical reaction network is one for which 
there exists a vector of positive species concentrations x* at which each of the reactions of the network is 
at equilibrium, that is, v{x*) = 0, see e.g. [29]. Such networks are necessarily reversible. Clearly every 
thermodynamically balanced network is complex-balanced. 

We now rewrite the dynamical equations for complex-balanced networks governed by mass action kinetics 
in terms of a known complex-equilibrium. It will be shown that such a form of equations has advantages in 
deriving stability properties of and also a model reduction method for complex-balanced networks. Recall 
equation ^ for general mass action reaction networks: 

x = -ZLExp (Z T Ln(x)) (8) 
Assume that the network is complex-balanced with a complex-equilibrium x* . Define 

K{x*) :=diagti(exp(ZfLn(x*))) 
where Zi denotes the i th column of Z. Equation Q can be rewritten as 

x = -ZLK(x*)K(x*)- 1 Exp(Z T Ln(x)) = -ZC(x*)Exp (z T Ln (^)) (9) 



where C(x*) := LK(x*). Note that since ijTL = 0, also l^£(x*) = 0. Furthermore since x* is a complex- 
equilibrium, we have 



£(x*)l c = £(x*)Exp(Z 7 Ln(— J J =0 

\ \x / / \ x=x * 

Hereafter, we refer to C(x*) as the weighted Laplacian of the graph of complexes associated with the given 
complex-balanced network. Both the row and column sums of the weighted Laplacian £{x*) are equal to 
zero It is this special property of the weighted Laplacian that we make use of in deriving all the results 
stated further on in this paper. 

3.3 The linkage classes of a graph of complexes 

A linkage class of a chemical reaction network is a maximal set of complexes {Ci, . . . ,Ck} such that Cj is 
connected by reactions to Cj for every i,jE{l,...,k},i^j. It can be easily verified that the number 
of linkage classes (I) of a network, which is equal to the number of connected components of the graph of 



2 In the literature on directed graphs (see e.g. [B]), C is called balanced. Note that the matrix L, having zero column sums 
but not zero row sums, is similar to the 'advection' set-up considered in j6]. 



complexes corresponding to the network, is given by £ = c— rank(i?) (the number of linkage classes in the 
terminology of \16\ [T2l I13j). The graph of complexes is connected, i.e., there is one linkage class in the 
network if and only if kex(B T ) = span(l c ). 

Assume that the reaction network has £ linkage classes. Assume that the i th linkage class has rj reactions 
between q complexes. Partition Z, B and C matrices according to the various linkage classes present in the 
network as follows: 



B 



Z\ Z2 

B x 

B 2 

... 

... 



.. Zi 




B t _! 









Br 



C(x*) = diag(£ 1 (x*),£ 2 (x*),...,C i - 1 (x*),C i (x*)) 



where for (i = !,...,£), Zi G M™ xc % B t G W lXTi and £i(x*) G R CiXCi denote the complex-stoichiometric 
matrix, incidence matrix and the weighted Laplacian matrices corresponding to equilibrium concentration 
x* respectively for the i th linkage class. Let Si denote the stoichiometric matrix of the i th linkage class. It 
is easy to see that Si = ZiB{. Observe that equation Q can be written as 



^Z i £ i (x*)Exp(zfLn(^ 



i=l 



(10) 



Remark 3.3. The law of conservation of mass states that there exists u G KT, such that Zfu G span (l Ci ) 
for i = 1, . . . ,£. This implies that u T x is a conserved moiety for the dynamics x = ZBv, for all forms of the 
reaction rate v(x). Indeed, Zju G span (l Ci ) implies u T ZB = 0, since Bft Ci = 0. 



4 Equilibria and stability of complex-balanced networks 

In this section, we make use of the compact mathematical formulation Q in order to derive properties of 
equilibria and stability of complex-balanced networks. 

4.1 Equilibria 

Our first result is a characterization of the set of all equilibria of a complex-balanced network in terms of a 
known equilibrium. 

Theorem 4.1. Consider a complex-balanced network governed by mass action kinetics. Let S G M mxr 
denote the stoichiometric matrix and assume that x* G WV: is a complex- equilibrium for the network. The 
following hold: 

1. x** G W± is another equilibrium for the network iff S T Ln (^r) = 0. 



2. Every equilibrium of the network is a complex- equilibrium. 

The proof of this theorem crucially makes use of the following lemma, which will also be the basis for 



the proof of Theorem 4.8 



Lemma 4.2. Let £{x*) be a balanced weighted Laplacian matrix as before. Then for any 7 € M c 

7 T £(x*)Exp( 7 ) > 
Moreover 

7 T £(a;*)Exp(7) = if and only if B T j = 

Proof. Let 7$ denote the z th element of 7 and let denote the negative of the element of C(x*) corresponding 
to the j th row and i th column. Note that kij > 0,i,j = 1, • • • , c. Using \^C{x*) = the expression 
— 7 T £(x*)Exp(7) can be rewritten as 

£i^i(7i-7l)fcii £#2(7i -12)1*21 ... £ i#c (7* - lc)k c i\ Exp( 7 ) 

c 

= - 7i)%iexp(7i) 

Furthermore, since the exponential function is strictly convex 

(/3 — a)exp(a) < exp(a) — exp(/3) 
for all ct,P, with equality if and only if a = /3. Hence 

c c 

- 7 T £(^*)Exp(7) = ^ J^(7i - 7 i )%exp(7 J -) < ^ ^ fyi(exp(7i) - exp( 7i )) 



'ID 



(12) 



£^ 1 fcii(exp(7i) - exp(7i)) 



1 



_£i^c ^(exp(7i) - exp( 7c )) 
-l^£( 3 ;*) T Exp(7) 
-(£(x*)l c ) T Exp( 7 ) = 0. 



since C{x*) is balanced. 



Furthermore, equality occurs in inequality ( 12 ) only when each of the terms within the summation on 



the left hand side is equal to the corresponding term within the summation on the right hand side. Since 
kij > 0, i 7^ j, if the i th complex reacts to the j th complex, it follows that 7^ = jj for each such i,j, which is 
equivalent to B T j = 0. ■ 



Proof, (of Theorem 4.1) The dynamics of the complex-balanced network with c complexes are given by 



x = -Z£(x*)Exp(Z T Ln 



where Z G M™ xc and C(x*) G M cxc are as defined in the previous section. Let B G M cxr denote the 
incidence matrix of the graph of complexes associated with the network. Assume that the reaction network 
has £ linkage classes. Assume that the i th linkage class has rj reactions between Cj complexes. Partition Z, 



B and C matrices according to the various linkage classes present in the network as in Section 3.3 Define 
Si := Z{Bi for % = 1, . . . , £. 

GO 

(Only If): Assume that x** is an equilibrium, that is 
-Z£(x*)Exp (Vim ( X ^)) = (13) 



Define 7 := Z T Ln (^*-)- Premultiplying equation (13) with Ln (^r), we get 
-7 T £(x*)Exp( 7 ) = 0. 



Hence by Lemma 4.2 B 7 = and thus 5 Ln f^nr) = 0. 



(//) Assume that S ,T Ln(^ r ) = 0. Hence for every linkage class i = 1, ...,£, Ln (^r) = 0, or, 
equivalently, Bfj = 0. This implies that ji = jj if the i th complex reacts to the j th complex or vice- versa. 
This in turn implies that for every linkage class i = 1, . . . ,£, Zj~Ln (^tt) consists of equal entries, or in other 
words it can be written as 

ZlLU (f^) = ^ 
where Tj G M. 

Since x* is a complex-equilibrium, C(x*)t c = 0. This implies that for i = 1, . . . ,£, Ci(x*)\ Ci = 0. Now, 



by evaluating the RHS of (10) at x**, we have 



Y^Zi&ix^Vxp (Vlu (^t)) = -^exp(r i )Zi£ i (x*)l Cj = 0. 



(<?.) Let x** denote an equilibrium as in the proof of the earlier part. Then S T Ln (^tt) = 0. We 
prove that x** is a complex-equilibrium. As shown earlier, Z?"Ln(^-) = Ljl^ for i = !,...,£. Since 
jCi (ac*)lc< = 0, this implies that (c.f., the discussion that preceeds Q and the form in (fTo))) 



-Bv(x**) = £(x*)Exp 



MS)) 



£i(x*)Exp(Z 1 T Ln(^)) 



£/(z*)Exp (Z £ T Ln(^)) 
From the above equation, it follows that x** is a complex-equilibrium. 



Remark 4.3. The steps followed in the proof of the above theorem are very similar to the proof of the 
characterization of the space of equilibria of a class of networks known as zero- deficiency networks presented 



in |13j . In the next subsection, we define zero-deficiency networks and prove that every zero-deficiency 



network that admits an equilibrium is complex-balanced. We emphasize here that the proof of Theorem 4.1 
that is presented in this paper is much more simple as compared to similar proofs provided in [J3] due to 
the use of the properties of the weighted Laplacian in the present manuscript. 

One may wonder to what extent the balanced weighted Laplacian matrix C(x*) depends on the choice 
of the complex-equilibrium x* . This dependency turns out to be very minor, strengthening the importance 
of this matrix for the analysis of the network. Indeed, consider any other complex-equilibrium x** . Then 
S T Lnx** = S T Lnx*, or equivalently B T Z T hnx** = B T Z T hnx* . Hence for the i-th connected component 
of the complex graph we have Bf ZfLnx** = Bf Z^hnx* , or equivalently, since ker Bf = spanl, 

Zfhnx** = ZfLnx* + q1 (14) 

for some constant Cj. Thus from the definition of £, it follows that Ci(x**) = diCi{x*) for some positive 
constant di. Hence, on every connected component of the graph of complexes, the balanced weighted 
Laplacian matrix C(x*) is unique up to multiplication by a positive constant. 

4.2 Zero-deficiency networks 

We now introduce the notion of zero-deficient chemical reaction networks mentioned in Remark 14.31 This 
notion was introduced in the work of Feinberg [11] in order to relate the stoichiometry of a given network 
to the structure of the associated graph of complexes. 

Definition 4.4. The deficiency 5 of a chemical reaction network with complex- stoichiometric matrix Z , 
incidence matrix B and stoichiometric matrix S is defined as 

5 := rank(S) - iank(ZB) = rank(5) - rank(S') > (15) 
A reaction network has zero- deficiency if 5 = 0. 

Note that zero-deficiency is equivalent with ker(Z)n im(B) = 0, or with the mapping 

Z:imBcl c ^r 
being injective. 



Remark 4.5. The deficiency of a chemical reaction network has been defined in a different way in |llj . 
Denote by I the number of linkage classes of a given chemical reaction network. Note that I = c — rank(l?) 



as explained in Section 3.3 In [TT], deficiency 5 is defined as 



S ■- c - £ - rank(S) (16) 



It is easy to see that definitions (15) and (16) are equivalent. 



We now prove that every zero-deficiency network that admits an equilibrium is complex-balanced. Con- 
sequently all the results that we state for a complex-balanced network also hold for a zero-deficiency network 
that admits an equilibrium. 

Lemma 4.6. // a chemical reaction network is zero-deficient and admits an equilibrium, then it is complex- 
balanced. 

Proof. Consider a zero-deficient network with complex-stoichiometric matrix Z G ]R™ XC and incidence matrix 
B G M cxr . Let x* G denote an equilibrium for the given network. Then Sv(x*) = ZBv(x*) = and 
hence by zero-deficiency Bv(x*) = 0. Consequently x* is a complex-equilibrium and the network is complex- 
balanced. ■ 

The above lemma has been stated and proved earlier in Theorem 4.1, p. 192] in a different and more 
lengthy manner. 

Remark 4.7. For McKeithan's network it is easily seen that Z itself is already injective, thus implying 
zero-deficiency. 

4.3 Asymptotic stability 

We now show global asymptotic stability of complex-balanced networks. 

Theorem 4.8. Consider a complex-balanced network with stoichiometric matrix S G M mxr , an equilibrium 
x* G and dynamics given by equation Q. Assume that the network obeys the law of conservation of 



mass stated in Remark 3.3. Then for every initial concentration vector xq G the species concentration 



vector x converges as t — )• oo to an element of the set 

8 := {x** G | S T Ln(x**) = S T Ln(x*)}. (17) 
Proof. Define 

G{x) = x T Ln (J^j + (x* - x) T l m (18) 

Observe that G(x*) = We prove that 

G{x)>0 Vx/x*, (19) 

and is proper, i.e., for every real c > the set {x G MIp | G(x) < c} is compact. Let X{ and x* denote the i th 
elements of x and x* respectively. From the strict concavity of the logarithmic function, 

z - ln(z) > 1 (20) 



Vz6 M+ with equality occuring only when z = 1. Putting z = in equation (20), we get 



- r, + .r. In ( -1 ] > 



G defined by ( 18 1 is a standard Lyapunov function used in chemical reaction network theory (see for example [13] 



This implies that 



G(x) = £ 



i=i 



x* - Xi + Xi In 



> 0. 



with equality occuring only when x = x* , thus proving (19). Properness of G is readily checked. 
We first prove that 
8G. 



G(x) := -z-(x)x < \/x G 
ox 



and 



G(x) = if and only if x G £ . 
Observe that 



(21) 



(22) 



G{x) 



-Ln 



x 



x 



Z T C(x*)Exp (Z T Ln 
we thus obtain 



Defining 7 := Z Ln (-^ 

= - 7 T £(j;*)Exp(7) 

and the statement follows from Lemma 14.21 

We now prove that M.™ is forward invariant with respect to ([9]). Assume by contradiction that this is 
not the case, and that X{(t) = for some t and i G {1, . . . , m}. Without loss of generality assume that 
the species with concentration x% is present in at least one complex which is involved in a reaction. Let Cj 
be the subset of complexes which contain Xj, and denote by Zc t the matrix containing the columns of Z 
corresponding to the complexes in C«. Hence all elements of the z-th row Z^ are different from zero. Then 
it follows that Y\T=i X ^ C = if C G Ci and thus 



±i(t) = -Z i L(x* 



.117 :•'•,(')' 



•jc c 



>0, 



where Z l is the i-th row vector of Z. This inequality is due to the fact that the terms corresponding to the 
positive i-th diagonal element of the weighted Laplacian matrix C(x*) are all zero, while there is at least 
one term corresponding to a non-zero, and therefore strictly negative, off-diagonal element of C(x*). It is 
easy to see that the inequality also holds at points arbitrarily close to the boundaries of the positive orthant 
Rip except the origin. Recall that according to the law of conservation of mass, there exists u G such 
that u T x is a conserved moiety. Consequently, starting from an initial concentration vector xo G M+, the 
state trajectory x(-) can not reach the origin. This implies that M.™ must be forward invariant with respect 
to @. 

Since G is proper (in M.™) and the state trajectory x(-) remains in M™, (21) implies that x(-) is bounded 
in M™. Therefore, boundedness of x(-), together with equations (21) and (22), imply that the species 



concentration x converges to an element of the set £ by an application of LaSalle's invariance principle. 



Remark 4.9. The crux of the proofs of Theorems 4.1 and 4.8 is the inequality 7 T £(x*)Exp(7) > 0, for all 7 G 



R. This inequality holds because of balancedness of C and we make use of the convexity of the exponential 
function in order to prove it. 



Remark 4.10. By proving that is forward invariant with respect to (|9j), we have actually proved the 
persistence conjecture stated in [U p. 1488] of complex-balanced networks obeying the law of conservation 
of mass. Persistence conjecture roughly corresponds to the requirement of nonextinction of species concen- 
tration of a complex-balanced network. Forward invariance of with respect to ^ has been proved in 
|27l Section VII] for the case when the complex-stoichiometric matrix Z is injective. Proving invariance of 
is an intricate problem for general mass action kinetics; see e.g. [4] and the references quoted therein. 

4.4 Equilibrium concentration corresponding to an initial concentration 

In this section, we show that corresponding to every positive stoichiometric compatibility class (see equa- 
tion Q for a definition) of a complex-balanced chemical reaction network, there exists a unique complex- 



equilibrium in £ defined by equation (17). The proof that we provide for this result is very similar to 
the proof of the zero-deficiency theorem provided in [13] and is based on the following proposition therein. 
Recall from the Introduction that the product x ■ z G M m is defined element- wise. 

Proposition 4.11. Let U be a linear subspace ofM. m , and let x*,xq G Then there is a unique element 
fi G U 1 - , such that (x* • Exp(/x) — xq) G U . 

Proof. See proof of [T31 Proposition B.l, pp. 361-363]. ■ 

Theorem 4.12. Consider the complex-balanced chemical reaction network with dynamics given by equation 
Q and equilibrium set £. Then for every xq G R™ there exists a unique x\ G £ nS Xo with S XQ given by 
and the solution trajectory x(-) of Q with initial condition x(0) = xq converges for t — > 00 to x\. 



Proof. With reference to Proposition |4 . 1 1 [ define U = imS, and observe that U ± = ker S T . By Proposition 



4.11 



there exists a unique [i G ker£ T such that x* ■ Exp(/i) — xq G imS. Define x± := x* ■ Exp(^) G M+. It 
follows that S T fi = S T Ln (%) = 0, i.e., x\ G £. Also, since x\ — xq G im S, x\ G S Xo which is an invariant set 



of the dynamics. Together with Theorem 4.8 it follows that the state trajectory x{-) with initial condition 



x(0) = xq converges to the equilibrium x\ G £. 



5 Model reduction 

For biochemical reaction networks, model-order reduction is still underdeveloped. The singular perturbation 
method and quasi steady-state approximation (QSSA) approach are the most commonly used techniques, 
where the reduced state contains a part of the metabolites of the full model. In the thesis by Hardin [14] . the 
QSSA approach is extended by considering higher-order approximation in the computation of quasi steady- 
state. Sunnaker et al. in [26] proposed a reduction method by identifying variables that can be lumped 



together and can be used to infer back the original state. In Prescott & Papachristodoulou [23J, a method to 
compute the upper-bound of the error is proposed based on sum-of-squares programming. The application 
of these techniques to general kinetics laws, such as Michaelis-Menten, poses a significant computational 
problem. 

In this section, we propose a novel and simple method for model reduction of complex-balanced chemical 
reaction networks governed by mass action kinetics. Our method is based on the Kron reduction method 
for model reduction of resistive electrical networks described in [18]; see also |28j . Moreover, the resulting 
reduced-order model retains the structure of the kinetics and gives result to a reduced biochemical reaction 
network, which enables a direct biochemical interpretation. 

5.1 Description of the method 

Our model reduction method is based on reduction of the graph of complexes associated with the network. 
It is based on the following result regarding Schur complements of weighted Laplacian matrices. 

Proposition 5.1. Consider a complex-balanced network with a complex- equilibrium x* G and weighted 
Laplacian matrix £(x*) G M cxc corresponding to the equilibrium x* . Let V denote the set of vertices of 
the graph of complexes associated with the network. Then for any subset of vertices V r C V , the Schur 
complement £{x*) of C{x*) with respect to the indices corresponding to V r satisfies the following properties: 

1. All diagonal elements of £{x*) are positive and off-diagonal elements are nonnegative. 

2. tTC(x*) = and C(x*)t £ = 0, where c := c - dim(V r ). 

Proof. (1.) Follows from the proof of [2U| Theorem 3.11]; see also [2H] for the case of symmetric £. 

(2.) Without loss of generality, assume that the last c — c rows and columns of C(x*) correspond to the 
vertex set V r . Consider a partition of C(x*) given by 



where £ n (x*) G M ex£ , C i2 {x*) G R 2x ( c ~ e ) , C 2 i(x*) G M^)^ and £ 2 2(>*) G r( c -c) x ( c -c>. By definition, 
C(x*) = £ n (x*) - £i 2 (x*)£ 22 (x*)- 1 £ 2 i(x*) 



C n (x*) C 12 (x*) 
C 2 i(x*) C 22 (x*) 



(23) 



Since tj£(x*) = 0, we obtain 

lf£n(x*) + lj_ £ £ 21 (x*) = 
tl£ 12 (x*) + l T c _zC 22 (x*) = 



Using the above equations, we get 




llXAiO*) - C 12 (x*)C 22 {x*)- l C 21 (x*)) 
-\ T c _,C 2l {x*) + lj_ £ £ 22 (x*)£ 22 (x*)- 1 £ 21 (x*) = 



In a similar way, it can be proved that C(x*)lc = 0. 



From the above result, it follows that C(x*) obeys all the properties of the weighted Laplacian of a 
complex-balanced chemical network corresponding to the graph of complexes with vertex set V — V r . Thus 



Proposition 5.1 can be directly applied to the graph of complexes, yielding a reduction of the chemical 
reaction network by reducing the number of complexes. Consider a complex-balanced reaction network 
described in the standard form ^ 



£ : 



-Z£(x*)Exp (Z T Ln 



Reduction will be performed by deleting certain complexes in the graph of complexes, resulting in a reduced 
graph of complexes with weighted Laplacian C(x*). Furthermore, leaving out the corresponding columns of 
the complex-stoichiometric matrix Z one obtains a reduced complex-stoichiometric matrix Z (with as many 
columns as the remaining number of complexes in the graph of complexes), leading to the reduced reaction 
network 



£ : 



-Z£(x*)Exp (z T Ln 



(24) 



Note that £ is again a complex-balanced chemical reaction network governed by mass action kinetics, with 
a reduced number of complexes and with, in general, a different set of reactions (edges of the graph of 
complexes). Furthermore, the complex-equilibrium x* of the original reaction network £ is a complex- 
equilibrium of the reduced network £ as well. 

An interpretation of the reduced network £ can be given as follows. Consider a subset V r of the set of 
all complexes, which we wish to leave out in the reduced network. Consider the partition of C(x*) as given 



by equation ( 23 ) and a corresponding partition of Z given by 
Z = Z~\ Zo 



(25) 



where V r corresponds to the last part of the indices (denoted by 2), in order to write out the dynamics of 
£ as 



x 



Z\ Z<i 



£n(x*) C l2 {x*) 
£21 (**) C 22 (x*) 



Exp(ZfLn(^)) 
Exp(Z 2 r Ln(^)) 



Consider now the auxiliary dynamical system 

C u (x*) C 12 (x*) 
C 2l (x*) C 22 (x*) 









w 2 




0. 



w 2 = -C 22 (x*) 1 C 2 i(x*)w 1 , 
leading to the reduced auxiliary dynamics defined by the Schur complement 



yi = -[Cn(x*) - C 12 {x*)C 22 (x*)- l £ 21 (x*)) Wl = -C(x*)wi 



reduced network E given in (24). 



Putting back in w\ = Exp ( Zfhn (^r) J , making use of x = Z\y\ + ^2^2 = Ziyi = Zyi, we then obtain the 



We derive the following properties relating £ and E. 
Proposition 5.2. Consider the complex-balanced reaction network £ and its reduced order model £ given 



by {24)- Denote their sets of equilibria by £, respectively £. Then £ C £. 



Proof. Let B denote the incidence matrix for E and let c and c denote the number of complexes in E and 
E respectively. Assume that the graph of complexes is connected; otherwise the same argument can be 
repeated for every component (linkage class). It follows that k.ev(B T ) = span(l c ). Let x** G £. We show 
that x** G £. Let Z and Z denote the complex-stoichiometric matrices of E and E respectively. Let C(x*) 
denote the weighted Laplacian of the graph of complexes of E corresponding to an equilibrium x* . Let C(x*) 
denote the Schur complement of C{x*) corresponding to the reduced model E. 

Since x** G £, B T Z T Ln (^) = 0. It follows that Z T Ln (^) G span(l c ). This implies that Z T Ln (^) G 
span(lg) since the columns of Z form a subset of the columns of Z. From Proposition 
that £(x*)Exp (Z T Ln (^-r) ) = 0. Consequently x** G £. This concludes the proof. 



5.1 



it now follows 



5.2 Effect of Model Reduction 

In this section, we show the effect of our model reduction method on two types of complex-balanced networks 
with a single linkage class. In other words, we give an interpretation of our reduced model in terms of its 
corresponding full model for two types of networks. Note that deletion of a set of complexes in one linkage 
class does not affect the reactions of the other linkage classes of the network. 

Type 1: 

fcl £,'2 ^3 kn-l 

Full Network: d ^ C 2 ^ C 3 ^ ^ C n (26) 

fc_l fc_2 fc_ 3 fc_( n _i) 



Reduced Network: C\ ^ C3 ^ f± C n (27) 

fc -i fc -2 fc-(n-l) 

These are complex-balanced networks with reversible reactions occuring between consecutive elements 



of the set of distinct complexes {Ci,C2, . . . ,C n } as in (26). The reduced network obtained by deleting the 



complex C2 can be computed to be (27). The two reactions, C\ ^ C2 and C2 f± C3 in the full network 
are replaced by one reaction C\ ^ C3 in the reduced network. This reaction is again a reversible reaction 



governed by mass action kinetics, with rate constants given by (27). 

The transient behaviour of the metabolites involved in the complexes of the reduced model will approxi- 
mately be the same as that of the full model if the metabolites involved in C2 reach their steady states much 
faster than the remaining metabolites. In this case, we can safely impose the condition that the metabolites 



involved in C2 are at constant concentration in order to obtain the reduced model (27) with similar transient 



behaviour as that of (26). The rule of induction may be applied in order to further reduce the model by 
deleting more complexes. 

A special case of Type 1 networks is C\ =± C2. Deletion of the complex C2 in this case is equivalent to 
deletion of the linkage class from the network. Such a deletion provides a close approximation to the original 
network if the reaction has very little effect on the dynamics of the network, i.e., if the reaction reaches its 
steady state much faster than the remaining reactions of the network. 




Figure 3: Type 2 full network Figure 4: Type 2 reduced network 



Type 2: These are complex-balanced networks between distinct complexes {C\,C2, ■ ■ ■ ,C n } as shown in 
Figure [3j McKeithan's network is an example of such networks. The reduced network obtained by deleting 
complex C2 is shown in Figure |4j Observe that the reaction C\ — > C$ in the reduced network has a different 
rate constant as compared to the full network. The three reactions C\ — >• C2, C2 — > C3 and C2 — > Co of the full 
network are replaced by one reaction C\ — > C3 in the reduced network. The rate constant for this reaction is 
given in Figure [4j All the remaining reactions of the reduced network occur in the same way as in the full 
network. 

In this case, the transient behaviour of the metabolites involved in the complexes of the reduced model 
will approximately be the same as that of the full model if the metabolites involved in C2 reach their steady 
states much faster than the remaining metabolites. Using the method described in this paper, we can study 
the effects of deleting other complexes like C\ or C n from the model. 

Observe that for all the types of networks discussed above, it is important to determine which of the 
complexes are to be deleted so that the reduced model approximates the full model reasonably well. 

Example 5.3. We have applied the model reduction method described in this paper to the model of T-cell 
interactions as in ©. We use the following numerical values: N = 19; 



kp,o = 


52; 


kp,i 


= 49; 


k p ,2 = 


41; 


k p ,3 = 


39; 


kp,4 — 


37; 


kp,5 


34; 


kpfi 


= 31; 


k p ,7 = 


29; 


kpfi = 


25; 


kp,g = 


19; 


■ = 


16; 


k p ,ii 


= 21; 


k p ,i2 = 


20; 


k P ,i3 = 


19; 


kp,14: = 


18; 




15; 


k P ,w 


= 24; 


k P ,n = 


13; 


k p ,i8 = 


= 7; 


k p ,i9 = 


= 5; 


k-1,0 = 


13; 


k-1,1 


= 29; 


k-1,2 = 


0.16; 


k-1,3 = 


1.4; 


fc-1,4 = 


2.3; 


k-1,5 = 


= 2; 


fe-1,6 : 


= 0.19; 


k-i,-j = 


0.33; 


k-1,8 = 


0.94; 


k-1,9 = 


0.67 


5-1,10 = 


0.31; 




= 0.21; 


fc-1,12 


= 3; 


fe-1,13 


= 5; 


fc-1,14 


= 1; 


fc-1,15 = 


= 11; 


fc-1,16 


= 0.8; 


fe-1,17 


= 7; 


fc-1,18 


= i; 


fc-1,19 = 


= 17. 
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Figure 5: Concentration profiles of T and M 



The initial value of each of the complexes Cj, i = 0, . . . , 19 is assumed to be equal to 0.01. The complexes 
T and M are assumed to have initial concentrations 1 and 2 respectively. We have performed model 
reduction by deleting 5 complexes C15, Cie, C17, Cig and C19. We have simulated the transient behaviour 
of the remaining complexes for the first two time units. The simulation results show that there is a good 
agreement between the transient behaviour of the concentration of most of such complexes when comparing 
the full network to the reduced network. Figure [5] depicts plots of comparison of the concentration profiles 
of T and M. 

An interpretation of model reduction of the above example model is as follows. By deleting complexes 
C15, Ci6, C17, Cis and C19, we assume that these complexes are at constant concentration. Since deleting 
these complexes results in a reduced model that closely mimics the original model, it follows that in the full 
model, these complexes reach an equilibrium much faster than the remaining complexes, so that assuming 
that these complexes are at constant concentrations results in a close approximation of the original model. 



6 Conclusion 

In this paper, we have provided a compact mathematical formulation for the dynamics of complex-balanced 
networks. We have made use of this formulation for the determination of equilibria and the asymptotic 
stability of such networks. The methods that have been employed are very similar to the ones used in |13j . 
but the difference is that our proofs are much more concise than the ones presented in [13J due to the use of 
properties of balanced weighted Laplacian matrices of complex-balanced networks. By proving the forward 
invariance of the positive orthant with respect to the dynamics of complex-balanced networks, we believe 
that we have proved the persistence conjecture for complex-balanced networks stated in [Ij. Furthermore, 
we have made use of the formulation in order to derive a model reduction technique for complex-balanced 
networks. 

A main challenge for further research is the extension of our results to chemical reaction networks with 
external fluxes and/or externally controlled concentrations. This will change the stability analysis consider- 
ably, due to the nonlinearity of the differential equations. Furthermore, it will also lead to scrutinizing the 



model reduction technique proposed in this paper from an external (input-output) point of view. 
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